<!DOCTYPE html>
<HTML lang = "en">
<HEAD>
  <meta charset="UTF-8"/>
  <meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
  <title>How Loops Work, An Introduction to Discrete Dynamics</title>
  

  <script type="text/x-mathjax-config">
    MathJax.Hub.Config({
      tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]},
      TeX: { equationNumbers: { autoNumber: "AMS" } }
    });
  </script>

  <script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
  </script>

  
<style>
pre.hljl {
    border: 1px solid #ccc;
    margin: 5px;
    padding: 5px;
    overflow-x: auto;
    color: rgb(68,68,68); background-color: rgb(251,251,251); }
pre.hljl > span.hljl-t { }
pre.hljl > span.hljl-w { }
pre.hljl > span.hljl-e { }
pre.hljl > span.hljl-eB { }
pre.hljl > span.hljl-o { }
pre.hljl > span.hljl-k { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kc { color: rgb(59,151,46); font-style: italic; }
pre.hljl > span.hljl-kd { color: rgb(214,102,97); font-style: italic; }
pre.hljl > span.hljl-kn { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kp { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kr { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kt { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-n { }
pre.hljl > span.hljl-na { }
pre.hljl > span.hljl-nb { }
pre.hljl > span.hljl-nbp { }
pre.hljl > span.hljl-nc { }
pre.hljl > span.hljl-ncB { }
pre.hljl > span.hljl-nd { color: rgb(214,102,97); }
pre.hljl > span.hljl-ne { }
pre.hljl > span.hljl-neB { }
pre.hljl > span.hljl-nf { color: rgb(66,102,213); }
pre.hljl > span.hljl-nfm { color: rgb(66,102,213); }
pre.hljl > span.hljl-np { }
pre.hljl > span.hljl-nl { }
pre.hljl > span.hljl-nn { }
pre.hljl > span.hljl-no { }
pre.hljl > span.hljl-nt { }
pre.hljl > span.hljl-nv { }
pre.hljl > span.hljl-nvc { }
pre.hljl > span.hljl-nvg { }
pre.hljl > span.hljl-nvi { }
pre.hljl > span.hljl-nvm { }
pre.hljl > span.hljl-l { }
pre.hljl > span.hljl-ld { color: rgb(148,91,176); font-style: italic; }
pre.hljl > span.hljl-s { color: rgb(201,61,57); }
pre.hljl > span.hljl-sa { color: rgb(201,61,57); }
pre.hljl > span.hljl-sb { color: rgb(201,61,57); }
pre.hljl > span.hljl-sc { color: rgb(201,61,57); }
pre.hljl > span.hljl-sd { color: rgb(201,61,57); }
pre.hljl > span.hljl-sdB { color: rgb(201,61,57); }
pre.hljl > span.hljl-sdC { color: rgb(201,61,57); }
pre.hljl > span.hljl-se { color: rgb(59,151,46); }
pre.hljl > span.hljl-sh { color: rgb(201,61,57); }
pre.hljl > span.hljl-si { }
pre.hljl > span.hljl-so { color: rgb(201,61,57); }
pre.hljl > span.hljl-sr { color: rgb(201,61,57); }
pre.hljl > span.hljl-ss { color: rgb(201,61,57); }
pre.hljl > span.hljl-ssB { color: rgb(201,61,57); }
pre.hljl > span.hljl-nB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nbB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nfB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nh { color: rgb(59,151,46); }
pre.hljl > span.hljl-ni { color: rgb(59,151,46); }
pre.hljl > span.hljl-nil { color: rgb(59,151,46); }
pre.hljl > span.hljl-noB { color: rgb(59,151,46); }
pre.hljl > span.hljl-oB { color: rgb(102,102,102); font-weight: bold; }
pre.hljl > span.hljl-ow { color: rgb(102,102,102); font-weight: bold; }
pre.hljl > span.hljl-p { }
pre.hljl > span.hljl-c { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-ch { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cm { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cp { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cpB { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cs { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-csB { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-g { }
pre.hljl > span.hljl-gd { }
pre.hljl > span.hljl-ge { }
pre.hljl > span.hljl-geB { }
pre.hljl > span.hljl-gh { }
pre.hljl > span.hljl-gi { }
pre.hljl > span.hljl-go { }
pre.hljl > span.hljl-gp { }
pre.hljl > span.hljl-gs { }
pre.hljl > span.hljl-gsB { }
pre.hljl > span.hljl-gt { }
</style>



  <style type="text/css">
  @font-face {
  font-style: normal;
  font-weight: 300;
}
@font-face {
  font-style: normal;
  font-weight: 400;
}
@font-face {
  font-style: normal;
  font-weight: 600;
}
html {
  font-family: sans-serif; /* 1 */
  -ms-text-size-adjust: 100%; /* 2 */
  -webkit-text-size-adjust: 100%; /* 2 */
}
body {
  margin: 0;
}
article,
aside,
details,
figcaption,
figure,
footer,
header,
hgroup,
main,
menu,
nav,
section,
summary {
  display: block;
}
audio,
canvas,
progress,
video {
  display: inline-block; /* 1 */
  vertical-align: baseline; /* 2 */
}
audio:not([controls]) {
  display: none;
  height: 0;
}
[hidden],
template {
  display: none;
}
a:active,
a:hover {
  outline: 0;
}
abbr[title] {
  border-bottom: 1px dotted;
}
b,
strong {
  font-weight: bold;
}
dfn {
  font-style: italic;
}
h1 {
  font-size: 2em;
  margin: 0.67em 0;
}
mark {
  background: #ff0;
  color: #000;
}
small {
  font-size: 80%;
}
sub,
sup {
  font-size: 75%;
  line-height: 0;
  position: relative;
  vertical-align: baseline;
}
sup {
  top: -0.5em;
}
sub {
  bottom: -0.25em;
}
img {
  border: 0;
}
svg:not(:root) {
  overflow: hidden;
}
figure {
  margin: 1em 40px;
}
hr {
  -moz-box-sizing: content-box;
  box-sizing: content-box;
  height: 0;
}
pre {
  overflow: auto;
}
code,
kbd,
pre,
samp {
  font-family: monospace, monospace;
  font-size: 1em;
}
button,
input,
optgroup,
select,
textarea {
  color: inherit; /* 1 */
  font: inherit; /* 2 */
  margin: 0; /* 3 */
}
button {
  overflow: visible;
}
button,
select {
  text-transform: none;
}
button,
html input[type="button"], /* 1 */
input[type="reset"],
input[type="submit"] {
  -webkit-appearance: button; /* 2 */
  cursor: pointer; /* 3 */
}
button[disabled],
html input[disabled] {
  cursor: default;
}
button::-moz-focus-inner,
input::-moz-focus-inner {
  border: 0;
  padding: 0;
}
input {
  line-height: normal;
}
input[type="checkbox"],
input[type="radio"] {
  box-sizing: border-box; /* 1 */
  padding: 0; /* 2 */
}
input[type="number"]::-webkit-inner-spin-button,
input[type="number"]::-webkit-outer-spin-button {
  height: auto;
}
input[type="search"] {
  -webkit-appearance: textfield; /* 1 */
  -moz-box-sizing: content-box;
  -webkit-box-sizing: content-box; /* 2 */
  box-sizing: content-box;
}
input[type="search"]::-webkit-search-cancel-button,
input[type="search"]::-webkit-search-decoration {
  -webkit-appearance: none;
}
fieldset {
  border: 1px solid #c0c0c0;
  margin: 0 2px;
  padding: 0.35em 0.625em 0.75em;
}
legend {
  border: 0; /* 1 */
  padding: 0; /* 2 */
}
textarea {
  overflow: auto;
}
optgroup {
  font-weight: bold;
}
table {
  font-family: monospace, monospace;
  font-size : 0.8em;
  border-collapse: collapse;
  border-spacing: 0;
}
td,
th {
  padding: 0;
}
thead th {
    border-bottom: 1px solid black;
    background-color: white;
}
tr:nth-child(odd){
  background-color: rgb(248,248,248);
}


/*
* Skeleton V2.0.4
* Copyright 2014, Dave Gamache
* www.getskeleton.com
* Free to use under the MIT license.
* http://www.opensource.org/licenses/mit-license.php
* 12/29/2014
*/
.container {
  position: relative;
  width: 100%;
  max-width: 960px;
  margin: 0 auto;
  padding: 0 20px;
  box-sizing: border-box; }
.column,
.columns {
  width: 100%;
  float: left;
  box-sizing: border-box; }
@media (min-width: 400px) {
  .container {
    width: 85%;
    padding: 0; }
}
@media (min-width: 550px) {
  .container {
    width: 80%; }
  .column,
  .columns {
    margin-left: 4%; }
  .column:first-child,
  .columns:first-child {
    margin-left: 0; }

  .one.column,
  .one.columns                    { width: 4.66666666667%; }
  .two.columns                    { width: 13.3333333333%; }
  .three.columns                  { width: 22%;            }
  .four.columns                   { width: 30.6666666667%; }
  .five.columns                   { width: 39.3333333333%; }
  .six.columns                    { width: 48%;            }
  .seven.columns                  { width: 56.6666666667%; }
  .eight.columns                  { width: 65.3333333333%; }
  .nine.columns                   { width: 74.0%;          }
  .ten.columns                    { width: 82.6666666667%; }
  .eleven.columns                 { width: 91.3333333333%; }
  .twelve.columns                 { width: 100%; margin-left: 0; }

  .one-third.column               { width: 30.6666666667%; }
  .two-thirds.column              { width: 65.3333333333%; }

  .one-half.column                { width: 48%; }

  /* Offsets */
  .offset-by-one.column,
  .offset-by-one.columns          { margin-left: 8.66666666667%; }
  .offset-by-two.column,
  .offset-by-two.columns          { margin-left: 17.3333333333%; }
  .offset-by-three.column,
  .offset-by-three.columns        { margin-left: 26%;            }
  .offset-by-four.column,
  .offset-by-four.columns         { margin-left: 34.6666666667%; }
  .offset-by-five.column,
  .offset-by-five.columns         { margin-left: 43.3333333333%; }
  .offset-by-six.column,
  .offset-by-six.columns          { margin-left: 52%;            }
  .offset-by-seven.column,
  .offset-by-seven.columns        { margin-left: 60.6666666667%; }
  .offset-by-eight.column,
  .offset-by-eight.columns        { margin-left: 69.3333333333%; }
  .offset-by-nine.column,
  .offset-by-nine.columns         { margin-left: 78.0%;          }
  .offset-by-ten.column,
  .offset-by-ten.columns          { margin-left: 86.6666666667%; }
  .offset-by-eleven.column,
  .offset-by-eleven.columns       { margin-left: 95.3333333333%; }

  .offset-by-one-third.column,
  .offset-by-one-third.columns    { margin-left: 34.6666666667%; }
  .offset-by-two-thirds.column,
  .offset-by-two-thirds.columns   { margin-left: 69.3333333333%; }

  .offset-by-one-half.column,
  .offset-by-one-half.columns     { margin-left: 52%; }

}
html {
  font-size: 62.5%; }
body {
  font-size: 1.5em; /* currently ems cause chrome bug misinterpreting rems on body element */
  line-height: 1.6;
  font-weight: 400;
  font-family: "Raleway", "HelveticaNeue", "Helvetica Neue", Helvetica, Arial, sans-serif;
  color: #222; }
h1, h2, h3, h4, h5, h6 {
  margin-top: 0;
  margin-bottom: 2rem;
  font-weight: 300; }
h1 { font-size: 3.6rem; line-height: 1.2;  letter-spacing: -.1rem;}
h2 { font-size: 3.4rem; line-height: 1.25; letter-spacing: -.1rem; }
h3 { font-size: 3.2rem; line-height: 1.3;  letter-spacing: -.1rem; }
h4 { font-size: 2.8rem; line-height: 1.35; letter-spacing: -.08rem; }
h5 { font-size: 2.4rem; line-height: 1.5;  letter-spacing: -.05rem; }
h6 { font-size: 1.5rem; line-height: 1.6;  letter-spacing: 0; }

p {
  margin-top: 0; }
a {
  color: #1EAEDB; }
a:hover {
  color: #0FA0CE; }
.button,
button,
input[type="submit"],
input[type="reset"],
input[type="button"] {
  display: inline-block;
  height: 38px;
  padding: 0 30px;
  color: #555;
  text-align: center;
  font-size: 11px;
  font-weight: 600;
  line-height: 38px;
  letter-spacing: .1rem;
  text-transform: uppercase;
  text-decoration: none;
  white-space: nowrap;
  background-color: transparent;
  border-radius: 4px;
  border: 1px solid #bbb;
  cursor: pointer;
  box-sizing: border-box; }
.button:hover,
button:hover,
input[type="submit"]:hover,
input[type="reset"]:hover,
input[type="button"]:hover,
.button:focus,
button:focus,
input[type="submit"]:focus,
input[type="reset"]:focus,
input[type="button"]:focus {
  color: #333;
  border-color: #888;
  outline: 0; }
.button.button-primary,
button.button-primary,
input[type="submit"].button-primary,
input[type="reset"].button-primary,
input[type="button"].button-primary {
  color: #FFF;
  background-color: #33C3F0;
  border-color: #33C3F0; }
.button.button-primary:hover,
button.button-primary:hover,
input[type="submit"].button-primary:hover,
input[type="reset"].button-primary:hover,
input[type="button"].button-primary:hover,
.button.button-primary:focus,
button.button-primary:focus,
input[type="submit"].button-primary:focus,
input[type="reset"].button-primary:focus,
input[type="button"].button-primary:focus {
  color: #FFF;
  background-color: #1EAEDB;
  border-color: #1EAEDB; }
input[type="email"],
input[type="number"],
input[type="search"],
input[type="text"],
input[type="tel"],
input[type="url"],
input[type="password"],
textarea,
select {
  height: 38px;
  padding: 6px 10px; /* The 6px vertically centers text on FF, ignored by Webkit */
  background-color: #fff;
  border: 1px solid #D1D1D1;
  border-radius: 4px;
  box-shadow: none;
  box-sizing: border-box; }
/* Removes awkward default styles on some inputs for iOS */
input[type="email"],
input[type="number"],
input[type="search"],
input[type="text"],
input[type="tel"],
input[type="url"],
input[type="password"],
textarea {
  -webkit-appearance: none;
     -moz-appearance: none;
          appearance: none; }
textarea {
  min-height: 65px;
  padding-top: 6px;
  padding-bottom: 6px; }
input[type="email"]:focus,
input[type="number"]:focus,
input[type="search"]:focus,
input[type="text"]:focus,
input[type="tel"]:focus,
input[type="url"]:focus,
input[type="password"]:focus,
textarea:focus,
select:focus {
  border: 1px solid #33C3F0;
  outline: 0; }
label,
legend {
  display: block;
  margin-bottom: .5rem;
  font-weight: 600; }
fieldset {
  padding: 0;
  border-width: 0; }
input[type="checkbox"],
input[type="radio"] {
  display: inline; }
label > .label-body {
  display: inline-block;
  margin-left: .5rem;
  font-weight: normal; }
ul {
  list-style: circle; }
ol {
  list-style: decimal; }
ul ul,
ul ol,
ol ol,
ol ul {
  margin: 1.5rem 0 1.5rem 3rem;
  font-size: 90%; }
li > p {margin : 0;}
th,
td {
  padding: 12px 15px;
  text-align: left;
  border-bottom: 1px solid #E1E1E1; }
th:first-child,
td:first-child {
  padding-left: 0; }
th:last-child,
td:last-child {
  padding-right: 0; }
button,
.button {
  margin-bottom: 1rem; }
input,
textarea,
select,
fieldset {
  margin-bottom: 1.5rem; }
pre,
blockquote,
dl,
figure,
table,
p,
ul,
ol,
form {
  margin-bottom: 1.0rem; }
.u-full-width {
  width: 100%;
  box-sizing: border-box; }
.u-max-full-width {
  max-width: 100%;
  box-sizing: border-box; }
.u-pull-right {
  float: right; }
.u-pull-left {
  float: left; }
hr {
  margin-top: 3rem;
  margin-bottom: 3.5rem;
  border-width: 0;
  border-top: 1px solid #E1E1E1; }
.container:after,
.row:after,
.u-cf {
  content: "";
  display: table;
  clear: both; }

pre {
  display: block;
  padding: 9.5px;
  margin: 0 0 10px;
  font-size: 13px;
  line-height: 1.42857143;
  word-break: break-all;
  word-wrap: break-word;
  border: 1px solid #ccc;
  border-radius: 4px;
}

pre.hljl {
  margin: 0 0 10px;
  display: block;
  background: #f5f5f5;
  border-radius: 4px;
  padding : 5px;
}

pre.output {
  background: #ffffff;
}

pre.code {
  background: #ffffff;
}

pre.julia-error {
  color : red
}

code,
kbd,
pre,
samp {
  font-family: Menlo, Monaco, Consolas, "Courier New", monospace;
  font-size: 0.9em;
}


@media (min-width: 400px) {}
@media (min-width: 550px) {}
@media (min-width: 750px) {}
@media (min-width: 1000px) {}
@media (min-width: 1200px) {}

h1.title {margin-top : 20px}
img {max-width : 100%}
div.title {text-align: center;}

  </style>
</HEAD>

<BODY>
  <div class ="container">
    <div class = "row">
      <div class = "col-md-12 twelve columns">
        <div class="title">
          <h1 class="title">How Loops Work, An Introduction to Discrete Dynamics</h1>
          <h5>Chris Rackauckas</h5>
          <h5>September 15th, 2020</h5>
        </div>

        <p>As we saw with the physics-informed neural networks, the basics of most scientific models are dynamical systems. Thus if we want to start to dig into deeper methods, we will need to start looking into the theory and practice of nonlinear dynamical systems. In this lecture we will go over the basic properties of dynamical systems and understand their general behavior through code. We will also learn the idea of stability as an asymtopic property of a mapping, and understand when a system is stable.</p>
<h2>Discrete Dynamical Systems</h2>
<p>A discrete dynamical system is a system which updates through discrete updates:</p>
<p class="math">\[
u_{n+1} = model(u_n,n)
\]</p>
<p>There are many examples of a discrete dynamical system found throughout the scientific literature. For example, many ecological models are discrete dynamical systems, with the most famous being the logistic map:</p>
<p class="math">\[
u_{n+1} = r u_n (1 - u_n)
\]</p>
<p>describing the growth of a population with a carrying capacity of 1 and a growth rate of <code>r</code>. Another way in which discrete dynamical systems are often encountered is through time series models. These are generally seen in financial forecasting and For example, the autoregressive model AR1 is the following linear dynamical system:</p>
<p class="math">\[
u_{n+1} = \alpha u_n + \epsilon_n
\]</p>
<p>where <span class="math">$\epsilon$</span> is a standard normal random number. The AR&#40;k&#41; model allows itself to update using delays as well:</p>
<p class="math">\[
u_{n+1} = \sum_{j=0}^{k-1} \alpha_j u_{n-j} + \epsilon_n
\]</p>
<p>The ARMA model is one that allows using delays on the randomness as well:</p>
<p class="math">\[
u_{n+1} = \sum_{j=0}^{k-1} (\alpha_j u_{n-j}  + \beta_j \epsilon_{n-j})
\]</p>
<p>Another embodiment of a discrete dynamical system is a Recurrent Neural Network &#40;RNN&#41;. In its simplest form, a RNN is a system of the form:</p>
<p class="math">\[
u_{n+1} = u_n + f(u_n,\theta)
\]</p>
<p>where <span class="math">$f$</span> is a neural network parameterized by <span class="math">$\theta$</span>.</p>
<p>Note that discrete dyamical systems are even more fundamental than just the ones shown. In any case where a continuous model is discretized to loop on the computer, the resulting algorithm is a discrete dynamical system and thus evolves according to its properties. This fact will be revisited later.</p>
<h2>Properties of Linear Dynamical Systems</h2>
<p>First let&#39;s take a look at the scalar linear dynamical system:</p>
<p class="math">\[
u_{n+1} = \alpha u_{n}
\]</p>
<p>We want to ask what the global or geometric behavior of this system is. We can do this by expanding out the system. Notice that if <span class="math">$u_0$</span> is known, then</p>
<p class="math">\[
u_{n+1} = \alpha^n u_0
\]</p>
<p>The global behavior can then be categorized as:</p>
<ul>
<li><p>If <span class="math">$\Vert \alpha \Vert < 1$</span>, then <span class="math">$u_n \rightarrow 0$</span></p>
</li>
<li><p>If <span class="math">$\Vert \alpha \Vert > 1$</span>, then <span class="math">$u_n \rightarrow \infty$</span></p>
</li>
</ul>
<p>If <span class="math">$\Vert \alpha \Vert = 1$</span>, then <span class="math">$u_n \rightarrow u_0$</span> if everything is in the real numbers, but more complex dynamics can occur on the complex plane.</p>
<h2>Nonlinear Geometric Dynamics</h2>
<p>The Geometric Theory of Dynamical Systems is the investigation of their long-term properties and the geometry of the phase space which they occupy. Let&#39;s start looking at this in practical terms: how do nonlinear update equations act as time goes to infinity?</p>
<h4>Banach Fixed Point Theorem</h4>
<p>There are surprisingly simple results that we can prove. First let&#39;s recall the Banach Fixed Point Theorem &#40;also known as the Contraction Mapping Theorem&#41;. Let <span class="math">$(X,d)$</span> be a metric space &#40;<span class="math">$X$</span> is the set of points we are thinking of, here the real numbers. <span class="math">$d$</span> is a distance function&#41;. <span class="math">$f$</span> is a contraction mapping if</p>
<p class="math">\[
d(f(x),f(y)) \leq q d(x,y)
\]</p>
<p>where <span class="math">$q < 1$</span>, that is, if applying <span class="math">$f$</span> always decreases the distance. The theorem then states that if <span class="math">$f$</span> is a contraction mapping, then there is a unique fixed point &#40;point <span class="math">$x^\ast$</span> where <span class="math">$f(x^\ast)=x^\ast$</span>&#41; and a sequence such that <span class="math">$x_0 \rightarrow x^\ast$</span> where</p>
<p class="math">\[
x_{n+1} = f(x_n)
\]</p>
<p>The proof is by induction, showing that the sequence is Cauchy. For some <span class="math">$m>n$</span> we do by the Triangle Inequality</p>
<p class="math">\[
d(x_m,x_n) \leq d(x_{m},x_{m-1}) + \ldots + d(x_{n+1},x_n)
\]</p>
<p>then apply the contraction relation down to the bottom:</p>
<p class="math">\[
d(x_m,x_n) \leq q^{m-1} d(x_{1},x_{0}) + \ldots + q^{n} d(x_{1},x_0)
\]</p>
<p class="math">\[
d(x_m,x_n) \leq q^n d(x_{1},x_{0}) \sum_{k=0}^{m-n-1} q^k
\]</p>
<p>and since adding more never hurts:</p>
<p class="math">\[
d(x_m,x_n) \leq q^n d(x_{1},x_{0}) \sum_{k=0}^{\infty} q^k
\]</p>
<p>But that summation is just a geometric series now, and since <span class="math">$q<1$</span> we know it converges to <span class="math">$1/(1-q)$</span>, and so we get:</p>
<p class="math">\[
d(x_m,x_n) \leq \frac{q^n}{1-q} d(x_{1},x_{0})
\]</p>
<p>The coefficient converges to zero as <span class="math">$n$</span> increases, and so the sequence must be Cauchy, which implies there&#39;s a unique fixed point.</p>
<h4>Stability of Linear Discrete Dynamical Systems</h4>
<p>Now let&#39;s take a mapping <span class="math">$f$</span> which is sufficiently nice &#40;<span class="math">$f \in C^1$</span>, i.e. the derivative of <span class="math">$f$</span> exists and is continuous&#41;, where</p>
<p class="math">\[
x_{n+1} = f(x_n)
\]</p>
<p>Assume that <span class="math">$\Vert f^\prime (x^\ast) \Vert < 1$</span> at some point where <span class="math">$f(x)=x$</span>. Then by continuity of the second derivative, it follows that there is a neighborhood where <span class="math">$\Vert f^\prime (x) \Vert < 1$</span> &#40;&#41;. Now recall that this means</p>
<p class="math">\[
\frac{df}{dx} \leq 1
\]</p>
<p>which means that, for any <span class="math">$x$</span> and <span class="math">$y$</span> in the neighborhood,</p>
<p class="math">\[
\Vert \frac{f(y)-f(x)}{y-x} \Vert \leq 1
\]</p>
<p>or</p>
<p class="math">\[
\Vert f(y)-f(x) \Vert \leq \Vert y-x \Vert
\]</p>
<p>This is essentially another way of saying that a function that is differentiable is Lipschitz, where we can use the derivative as the Lipschitz bound. But notice this means that, in this neighborhood, a function with a derivative less than 1 is a contraction mapping, and thus there is a limiting sequence which goes to the fixed point by the Banach Fixed Point Theorem. Furthermore, the uniquess guerentees that there is only one fixed point in a sufficiently small neighborhood where the derivative is all less than 1.</p>
<p>A way to interpret this result is that, any nice enough function <span class="math">$f$</span> is locally linear. Thus we can understand the global properties of <span class="math">$f$</span> by looking the linearization of its dynamics, where the best linear approximation is the linear function <span class="math">$f^\prime (x) x$</span>. This means that we can think of</p>
<p class="math">\[
x_{n+1} = f(x_n)
\]</p>
<p>locally as being approximated by</p>
<p class="math">\[
x_{n+1} = f^\prime (x) x_n
\]</p>
<p>and so if the derivative is less than 1 in some neighborhood of a fixed point, then we locally have a linear dynamical system which looks like the simple <span class="math">$x_{n+1} = \alpha x_n$</span> where <span class="math">$\alpha <1$</span>, and so we get the same convergence property.</p>
<p>This is termed &quot;stability&quot; since, if you are a little bit off from the fixed point, you will tend to go right back to it. An unstable fixed point is one where you fall away. And what happens when the derivative is one? There are various forms of semi-stability that can be proved which go beyond the topic of this course.</p>
<h4>Update Form</h4>
<p>Now let&#39;s look at another form:</p>
<p class="math">\[
x_{n+1} = x_n + f(x_n)
\]</p>
<p>For example, this is what we generally see with the recurrent neural network &#40;or, as we will find out later, this is how discretizations of continuous systems tend to look&#33;&#41;. In this case, we can say that this is a dynamical system</p>
<p class="math">\[
x_{n+1} = g(x_n)
\]</p>
<p>and so if <span class="math">$-2 < f^\prime < 0$</span>, then <span class="math">$\Vert g^\prime \Vert = \Vert 1 + f^\prime \Vert < 1$</span> and so we have the same stability idea except now with a condition shifted to zero instead of one.</p>
<h2>Multivariable Systems</h2>
<p>Now let <span class="math">$x \in R^k$</span> be a vector, and define discrete mappings:</p>
<p class="math">\[
x_{n+1} = f(x_n)
\]</p>
<p>To visualize this, let&#39;s write out the version for <span class="math">$x \in R^3$</span>:</p>
<p class="math">\[
x_{n+1}=\left[\begin{array}{c}
a_{n+1}\\
b_{n+1}\\
c_{n+1}
\end{array}\right]=\left[\begin{array}{c}
f_{1}(a_{n},b_{n},c_{n})\\
f_{2}(a_{n},b_{n},c_{n})\\
f_{3}(a_{n},b_{n},c_{n})
\end{array}\right]=f(x_{n})
\]</p>
<p>The linear multidimensional discrete dynamical system is:</p>
<p class="math">\[
x_{n+1} = A x_n
\]</p>
<p>The easiest way to analyze a multidimensional system is to turn it into a bunch of single dimension systems. To do this, assume that <span class="math">$A$</span> is diagonal. This means that there exists a diagonalization <span class="math">$A =P^{-1}DP$</span> where <span class="math">$P$</span> is the matrix of eigenvectors and <span class="math">$D$</span> is the diagonal matrix of eigenvalues. We can then decompose the system as follows:</p>
<p class="math">\[
Px_{n+1} = DPx_n
\]</p>
<p>and now define new variables <span class="math">$z_n = Px_n$</span>. In these variables,</p>
<p class="math">\[
z_{n+1} = D z_n
\]</p>
<p>but <span class="math">$D$</span> is diagonal, so this is a system of <span class="math">$k$</span> independent linear dynamical systems. We know that the linear dynamical system will converge to zero if <span class="math">$\Vert D_i \Vert < 1$</span>, and so this means that <span class="math">$z_n$</span> converges to zero if all of the eigenvalues are within the unit circle. Since <span class="math">$P0 = 0$</span>, this implies that if all of the eigenvalues of <span class="math">$A$</span> are in the unit circle, then <span class="math">$x_n \rightarrow 0$</span>.</p>
<p>A multidimensional version of the contraction mapping theorem is then proven exactly in this manner, meaning that if <span class="math">$f(x) = x$</span> and all eigenvalues of the Jacobian matrix &#40;the linearization of <span class="math">$f$</span>&#41; are in the unit circle, then <span class="math">$x$</span> is a unique fixed point in some neighborhood.</p>
<h4>Understanding Delayed Systems</h4>
<p>A similar property holds in linear dynamical systems with delays. Take</p>
<p class="math">\[
x_{n+1} = \sum_{j=0}^{k-1} \alpha_j x_{n-j}
\]</p>
<p>Notice that we can write this as a multidimensional non-delayed system. Let <span class="math">$x_n^i$</span> be the <span class="math">$i$</span>th term in the vector of the <span class="math">$n$</span> time. Then we have:</p>
<p class="math">\[
x_{n+1}^1 = \sum_{j=1}^{k-1} \alpha_{j-1} x_{n}^{j}
\]</p>
<p>as an equivalent way to write this, where</p>
<p class="math">\[
x_{n+1}^j = x_n^{j-1}
\]</p>
<p>for all of the other terms. Essentially, instead of a system with a delay, we store the memory in other terms of the vector, and keep shifting them down. However, this makes our system much easier to analyze. Instead of a linear delayed dynamical system, this is now a linear multidimensional dynamical system. Its characteristic polynomial is</p>
<p class="math">\[
\varphi(x) = 1 - \sum_{j=0}^{k-1} \alpha_j x^j
\]</p>
<p>and so if all of the roots are in the unit circle then this system is stable.</p>
<h4>Stochastic Dynamical Systems</h4>
<p>Now let&#39;s take a look again at the autoregressive process from time series analysis:</p>
<p class="math">\[
u_{n+1} = \sum_{j=0}^{k-1} \alpha_j u_{n-j} + \epsilon_n
\]</p>
<p>In a very quick handwavy way, we can understand such a system by seeing how the perturbations propogate. If <span class="math">$u_0 = 0$</span>, then the starting is just <span class="math">$\epsilon_0$</span>. If we assume all other <span class="math">$\epsilon_i = 0$</span>, then this system is the same as a linear dynamical system with delays. If all of the roots are in the unit circle, then it goes to zero, meaning the perturbation is forgotten or squashed over time.</p>
<p>We can analyze this more by using the moments. Notice that, by the linearity of the expected value,</p>
<p class="math">\[
\mathbb{E}[u_{n+1}] = \sum_{j=0}^{k-1} \alpha_j \mathbb{E}[u_{n-j}]
\]</p>
<p>is a deterministic linear dynamical system which converges if the roots are in the unit circle. This means that the mean stabilizes over time if all of the roots are in the unit circle. In time series analysis, this is called stationarity of the time series.</p>
<p>We can then also look at the stability of the variance as well. Recall that</p>
<p class="math">\[
\mathbb{V}[x] = \mathbb{E}[x^2] - \mathbb{E}[x]^2
\]</p>
<p>and so therefore</p>
<p class="math">\[
\mathbb{E}[u_{n+1}^2] = \mathbb{E}[\sum_{j=0}^{k-1} \alpha_j u_{n-j}^2]
\]</p>
<p>and with a bunch of analysis here, working in the same way with the same basic ideas, we can determine conditions on which the variance goes to zero.</p>
<h2>Periodicity and Chaos</h2>
<p>Stability is the simplest geometric dynamical property, but there are many others. For example, maps can also have periodic orbits, like:</p>
<p class="math">\[
u_{n+1} = -u_n
\]</p>
<p>will bounce back and forth between two values. These periodic orbits themselves have geometric properties, such as whether it&#39;s a stable periodic orbit &#40;points nearby are attracted to the periodic orbit&#41;. Periodic orbits have a length as well: this was a periodic orbit of length 2.</p>
<p>Chaos is another interesting property of a discrete dynamical system. It can be interpreted as a periodic orbit where the length is infinity. This can happen if, by changing a parameter, a period 2 orbit becomes a period 4, then a period 8, etc. &#40;a phonomenon known as period doubling&#41;, and when it goes beyond the accumulation point the &quot;infinite period orbit&quot; is reached and chaos is found. A homework problem will delve into the properties of chaos as an example of a simple embaressingly data-parallel problem.</p>
<h2>Efficient Implmentation of Dynamical Systems</h2>
<p>Dynamical systems are just loops, so the implementation is easy to understand. However, there are a few things that one must keep in mind in order to allow for efficient implementations.</p>
<h4>Higher order functions</h4>
<p>Functions which compute the solutions to dynamical systems are inherently <em>higher order functions</em>, which means it&#39;s a function which takes in a function as an argument. The following is a quick implementation:</p>


<pre class='hljl'>
<span class='hljl-s'>&quot;&quot;&quot;
`solve_system(f,u0,n)`

Solves the dynamical system

``u_{n+1} = f(u_n)``

for N steps. Returns the solution at step `n` with parameters `p`.

&quot;&quot;&quot;</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u0</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-n'>n</span><span class='hljl-oB'>-</span><span class='hljl-ni'>1</span><span class='hljl-t'>
    </span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span>
</pre>


<pre class="output">
Main.##WeaveSandBox#2406.solve_system
</pre>


<p>Notice the <code>&quot;&quot;&quot;</code> before the function: this is a docstring. When the Julia REPL is queried with <code>?solve_system</code> this will be the description that is displayed.</p>
<p>Now, is this function going to be efficient? Recall from the earlier discussion that:</p>
<ul>
<li><p>Type-stability is necessary for inference to carry forward type information.</p>
</li>
<li><p>Julia auto-specializes on input types.</p>
</li>
<li><p>Inlining can occur automatically for sufficiently small functions.</p>
</li>
</ul>
<p>From this information, we know that in order for this to be efficient, we require that the type of <code>f&#40;u&#41;</code> is inferred. But if <code>f</code> is a variable, how can that be inferred? In order for that to occur, we would have to know what <code>f</code> is, since not all functions will give the same output type. Additionally, in order to inline the function, we will have to know what the function is at compile-time. So, is it possible to make this implementation efficient?</p>
<p>It turns out that this does optimize due to one fact: every function is given by its own type. We can verify this by defining a function and checking its type:</p>


<pre class='hljl'>
<span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>p</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-t'>
</span><span class='hljl-nf'>typeof</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
typeof&#40;Main.##WeaveSandBox#2406.f&#41;
</pre>


<p>It displays <code>typeof&#40;f&#41;</code> to indicate that the function <code>f</code> is its own type, and thus at automatic specialization time the compiler knows all of the information about the function and thus inlines and performs inference correctly.</p>
<p>Note that this does mean that the function will need to recompile for every new <code>f</code>. This is similar to statically compiling a function for use in a C/Fortran library. What is the equivalent to using a function like a shared library or a shared object? This is given by FunctionWrappers.jl. This directly stores the function pointer in an object that can have shared type information in order to keep every function as the same type. However, the wrapped function has more information about the pointer... what&#39;s necessary?</p>
<p>The answer is that FunctionWrappers.jl allows for specifying the input and output types, in order for the wrapper to do the right assertions for inference to carry forward type stability, since in this case inference is not able to step through the function pointer.</p>
<h4>Quick Check</h4>
<p>What will approximately be the value of this dynamical system after 1000 steps if you start at <code>1.0</code> with parameter <code>p&#61;0.25</code>? Can you guess without solving the system? Think about steady states and stabiltiy.</p>


<pre class='hljl'>
<span class='hljl-nf'>solve_system</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.25</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
0.0
</pre>


<p>The answer is that it goes to zero. The steady states are the zeros of the polynomial, which are <code>0</code> and <code>p&#43;1</code>. It&#39;s reasonable to believe that it either goes to one of those 2 values or infinity. In the first step, <code>1^2 - 0.25 &#61; 0.75 &lt; 1</code> which suggests &#40;but doesn&#39;t confirm&#33;&#41; that it&#39;s a contraction. Notice that the derivative is <code>2u-p</code>, and so <code>u&#61;p&#43;1&#61;1.25</code> is not a stable steady state, and thus we go to zero. In fact, we can check a few values:</p>


<pre class='hljl'>
<span class='hljl-nf'>solve_system</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1.1</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.25</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
0.0
</pre>



<pre class='hljl'>
<span class='hljl-nf'>solve_system</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1.22</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.25</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
0.0
</pre>



<pre class='hljl'>
<span class='hljl-nf'>solve_system</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1.25</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.25</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
1.25
</pre>



<pre class='hljl'>
<span class='hljl-nf'>solve_system</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1.251</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.25</span><span class='hljl-p'>,</span><span class='hljl-ni'>20</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
Inf
</pre>


<p>Notice that the moment we go above the steady state <code>p&#43;1</code>, we exponentially grow to infinity.</p>
<p>Just to double check the implementation:</p>


<pre class='hljl'>
<span class='hljl-nf'>solve_system</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1.251</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.25</span><span class='hljl-p'>,</span><span class='hljl-ni'>10</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nf'>solve_system</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1.251</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.25</span><span class='hljl-p'>,</span><span class='hljl-ni'>100</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nf'>solve_system</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1.251</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.25</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
NaN
</pre>


<p>Those allocations are just the output, and notice it&#39;s independent of the loop count.</p>
<h4>Multidimensional System Implementations</h4>
<p>When we go to multidimensional systems, some care needs to be taken to decrease the number of allocations which are occuring. One of the ways to do this is to utilize statically sized arrays. For example, let&#39;s look at a discretization of the Lorenz system:</p>


<pre class='hljl'>
<span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>lorenz</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>α</span><span class='hljl-p'>,</span><span class='hljl-n'>σ</span><span class='hljl-p'>,</span><span class='hljl-n'>ρ</span><span class='hljl-p'>,</span><span class='hljl-n'>β</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>p</span><span class='hljl-t'>
  </span><span class='hljl-n'>du1</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-n'>α</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>σ</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-oB'>-</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]))</span><span class='hljl-t'>
  </span><span class='hljl-n'>du2</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-n'>α</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>ρ</span><span class='hljl-oB'>-</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>])</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>])</span><span class='hljl-t'>
  </span><span class='hljl-n'>du3</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-n'>α</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>β</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>])</span><span class='hljl-t'>
  </span><span class='hljl-p'>[</span><span class='hljl-n'>du1</span><span class='hljl-p'>,</span><span class='hljl-n'>du2</span><span class='hljl-p'>,</span><span class='hljl-n'>du3</span><span class='hljl-p'>]</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-n'>p</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-nfB'>0.02</span><span class='hljl-p'>,</span><span class='hljl-nfB'>10.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>28.0</span><span class='hljl-p'>,</span><span class='hljl-ni'>8</span><span class='hljl-oB'>/</span><span class='hljl-ni'>3</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nf'>solve_system</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
3-element Array&#123;Float64,1&#125;:
  1.4744010677851374
  0.8530017039412324
 20.62004063423844
</pre>


<p>Let&#39;s see what this gives us by saving:</p>


<pre class='hljl'>
<span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Vector</span><span class='hljl-p'>{</span><span class='hljl-nf'>typeof</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>)}(</span><span class='hljl-n'>undef</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u0</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-n'>n</span><span class='hljl-oB'>-</span><span class='hljl-ni'>1</span><span class='hljl-t'>
    </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-oB'>+</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-n'>to_plot</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
1000-element Array&#123;Array&#123;Float64,1&#125;,1&#125;:
 &#91;1.0, 0.0, 0.0&#93;
 &#91;0.8, 0.56, 0.0&#93;
 &#91;0.752, 0.9968000000000001, 0.008960000000000001&#93;
 &#91;0.80096, 1.3978492416000001, 0.023474005333333336&#93;
 &#91;0.92033784832, 1.8180538219817644, 0.04461448495326095&#93;
 &#91;1.099881043052353, 2.296260732619613, 0.07569952060880669&#93;
 &#91;1.339156980965805, 2.864603692722823, 0.12217448583728006&#93;
 &#91;1.6442463233172087, 3.5539673118971193, 0.19238159391549564&#93;
 &#91;2.026190521033191, 4.397339452147425, 0.2989931959555302&#93;
 &#91;2.5004203072560376, 5.431943011293093, 0.4612438424853632&#93;
 ⋮
 &#91;6.8089180814322185, 0.8987564841782779, 31.6759436385101&#93;
 &#91;5.6268857619814305, 0.3801973723631693, 30.108951163308078&#93;
 &#91;4.577548084057778, 0.13525687944525802, 28.545926978224173&#93;
 &#91;3.6890898431352737, 0.08257160199224252, 27.035860436772758&#93;
 &#91;2.9677861949066675, 0.15205611935372762, 25.600040161309696&#93;
 &#91;2.4046401797960795, 0.2914663505185634, 24.24373008707723&#93;
 &#91;1.9820054139405763, 0.46628657468365653, 22.964748583050085&#93;
 &#91;1.6788616460891923, 0.6565587545689172, 21.758445642263496&#93;
 &#91;1.4744010677851374, 0.8530017039412324, 20.62004063423844&#93;
</pre>



<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>Plots</span><span class='hljl-t'>
</span><span class='hljl-n'>x</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-n'>to_plot</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>][</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-nf'>length</span><span class='hljl-p'>(</span><span class='hljl-n'>to_plot</span><span class='hljl-p'>)]</span><span class='hljl-t'>
</span><span class='hljl-n'>y</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-n'>to_plot</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>][</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-nf'>length</span><span class='hljl-p'>(</span><span class='hljl-n'>to_plot</span><span class='hljl-p'>)]</span><span class='hljl-t'>
</span><span class='hljl-n'>z</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-n'>to_plot</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>][</span><span class='hljl-ni'>3</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-nf'>length</span><span class='hljl-p'>(</span><span class='hljl-n'>to_plot</span><span class='hljl-p'>)]</span><span class='hljl-t'>
</span><span class='hljl-nf'>plot</span><span class='hljl-p'>(</span><span class='hljl-n'>x</span><span class='hljl-p'>,</span><span class='hljl-n'>y</span><span class='hljl-p'>,</span><span class='hljl-n'>z</span><span class='hljl-p'>)</span>
</pre>


<img src=""  />

<p>This is the chaotic Lorenz attractor plotted in phase space, i.e. the values of the variables against each other.</p>
<p>Let&#39;s look at the implementation a little bit more. <code>u &#61; Vector&#123;typeof&#40;u0&#41;&#125;&#40;undef,n&#41;</code> is type-generic, meaning any <code>u0</code> can be used with that code. However, as a vector of vectors, it is a vector of pointers to contiguous memory, instead of being contiguous itself. Note that that means there is not much of a cost by not pre-specifying the size up front:</p>


<pre class='hljl'>
<span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save_push</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Vector</span><span class='hljl-p'>{</span><span class='hljl-nf'>typeof</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>)}(</span><span class='hljl-n'>undef</span><span class='hljl-p'>,</span><span class='hljl-ni'>1</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u0</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-n'>n</span><span class='hljl-oB'>-</span><span class='hljl-ni'>1</span><span class='hljl-t'>
    </span><span class='hljl-nf'>push!</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>))</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-nd'>@time</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
0.000062 seconds &#40;1.00 k allocations: 117.312 KiB&#41;
1000-element Array&#123;Array&#123;Float64,1&#125;,1&#125;:
 &#91;1.0, 0.0, 0.0&#93;
 &#91;0.8, 0.56, 0.0&#93;
 &#91;0.752, 0.9968000000000001, 0.008960000000000001&#93;
 &#91;0.80096, 1.3978492416000001, 0.023474005333333336&#93;
 &#91;0.92033784832, 1.8180538219817644, 0.04461448495326095&#93;
 &#91;1.099881043052353, 2.296260732619613, 0.07569952060880669&#93;
 &#91;1.339156980965805, 2.864603692722823, 0.12217448583728006&#93;
 &#91;1.6442463233172087, 3.5539673118971193, 0.19238159391549564&#93;
 &#91;2.026190521033191, 4.397339452147425, 0.2989931959555302&#93;
 &#91;2.5004203072560376, 5.431943011293093, 0.4612438424853632&#93;
 ⋮
 &#91;6.8089180814322185, 0.8987564841782779, 31.6759436385101&#93;
 &#91;5.6268857619814305, 0.3801973723631693, 30.108951163308078&#93;
 &#91;4.577548084057778, 0.13525687944525802, 28.545926978224173&#93;
 &#91;3.6890898431352737, 0.08257160199224252, 27.035860436772758&#93;
 &#91;2.9677861949066675, 0.15205611935372762, 25.600040161309696&#93;
 &#91;2.4046401797960795, 0.2914663505185634, 24.24373008707723&#93;
 &#91;1.9820054139405763, 0.46628657468365653, 22.964748583050085&#93;
 &#91;1.6788616460891923, 0.6565587545689172, 21.758445642263496&#93;
 &#91;1.4744010677851374, 0.8530017039412324, 20.62004063423844&#93;
</pre>



<pre class='hljl'>
<span class='hljl-nd'>@time</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save_push</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
0.025693 seconds &#40;19.79 k allocations: 1.150 MiB&#41;
1000-element Array&#123;Array&#123;Float64,1&#125;,1&#125;:
 &#91;1.0, 0.0, 0.0&#93;
 &#91;0.8, 0.56, 0.0&#93;
 &#91;0.752, 0.9968000000000001, 0.008960000000000001&#93;
 &#91;0.80096, 1.3978492416000001, 0.023474005333333336&#93;
 &#91;0.92033784832, 1.8180538219817644, 0.04461448495326095&#93;
 &#91;1.099881043052353, 2.296260732619613, 0.07569952060880669&#93;
 &#91;1.339156980965805, 2.864603692722823, 0.12217448583728006&#93;
 &#91;1.6442463233172087, 3.5539673118971193, 0.19238159391549564&#93;
 &#91;2.026190521033191, 4.397339452147425, 0.2989931959555302&#93;
 &#91;2.5004203072560376, 5.431943011293093, 0.4612438424853632&#93;
 ⋮
 &#91;6.8089180814322185, 0.8987564841782779, 31.6759436385101&#93;
 &#91;5.6268857619814305, 0.3801973723631693, 30.108951163308078&#93;
 &#91;4.577548084057778, 0.13525687944525802, 28.545926978224173&#93;
 &#91;3.6890898431352737, 0.08257160199224252, 27.035860436772758&#93;
 &#91;2.9677861949066675, 0.15205611935372762, 25.600040161309696&#93;
 &#91;2.4046401797960795, 0.2914663505185634, 24.24373008707723&#93;
 &#91;1.9820054139405763, 0.46628657468365653, 22.964748583050085&#93;
 &#91;1.6788616460891923, 0.6565587545689172, 21.758445642263496&#93;
 &#91;1.4744010677851374, 0.8530017039412324, 20.62004063423844&#93;
</pre>


<p>The first time Julia compiles the function, and the second is a straight call.</p>


<pre class='hljl'>
<span class='hljl-nd'>@time</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save_push</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
0.000112 seconds &#40;1.01 k allocations: 125.781 KiB&#41;
1000-element Array&#123;Array&#123;Float64,1&#125;,1&#125;:
 &#91;1.0, 0.0, 0.0&#93;
 &#91;0.8, 0.56, 0.0&#93;
 &#91;0.752, 0.9968000000000001, 0.008960000000000001&#93;
 &#91;0.80096, 1.3978492416000001, 0.023474005333333336&#93;
 &#91;0.92033784832, 1.8180538219817644, 0.04461448495326095&#93;
 &#91;1.099881043052353, 2.296260732619613, 0.07569952060880669&#93;
 &#91;1.339156980965805, 2.864603692722823, 0.12217448583728006&#93;
 &#91;1.6442463233172087, 3.5539673118971193, 0.19238159391549564&#93;
 &#91;2.026190521033191, 4.397339452147425, 0.2989931959555302&#93;
 &#91;2.5004203072560376, 5.431943011293093, 0.4612438424853632&#93;
 ⋮
 &#91;6.8089180814322185, 0.8987564841782779, 31.6759436385101&#93;
 &#91;5.6268857619814305, 0.3801973723631693, 30.108951163308078&#93;
 &#91;4.577548084057778, 0.13525687944525802, 28.545926978224173&#93;
 &#91;3.6890898431352737, 0.08257160199224252, 27.035860436772758&#93;
 &#91;2.9677861949066675, 0.15205611935372762, 25.600040161309696&#93;
 &#91;2.4046401797960795, 0.2914663505185634, 24.24373008707723&#93;
 &#91;1.9820054139405763, 0.46628657468365653, 22.964748583050085&#93;
 &#91;1.6788616460891923, 0.6565587545689172, 21.758445642263496&#93;
 &#91;1.4744010677851374, 0.8530017039412324, 20.62004063423844&#93;
</pre>


<p>or we can use <code>@btime</code>:</p>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>BenchmarkTools</span><span class='hljl-t'>
</span><span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
38.700 μs &#40;1001 allocations: 117.31 KiB&#41;
1000-element Array&#123;Array&#123;Float64,1&#125;,1&#125;:
 &#91;1.0, 0.0, 0.0&#93;
 &#91;0.8, 0.56, 0.0&#93;
 &#91;0.752, 0.9968000000000001, 0.008960000000000001&#93;
 &#91;0.80096, 1.3978492416000001, 0.023474005333333336&#93;
 &#91;0.92033784832, 1.8180538219817644, 0.04461448495326095&#93;
 &#91;1.099881043052353, 2.296260732619613, 0.07569952060880669&#93;
 &#91;1.339156980965805, 2.864603692722823, 0.12217448583728006&#93;
 &#91;1.6442463233172087, 3.5539673118971193, 0.19238159391549564&#93;
 &#91;2.026190521033191, 4.397339452147425, 0.2989931959555302&#93;
 &#91;2.5004203072560376, 5.431943011293093, 0.4612438424853632&#93;
 ⋮
 &#91;6.8089180814322185, 0.8987564841782779, 31.6759436385101&#93;
 &#91;5.6268857619814305, 0.3801973723631693, 30.108951163308078&#93;
 &#91;4.577548084057778, 0.13525687944525802, 28.545926978224173&#93;
 &#91;3.6890898431352737, 0.08257160199224252, 27.035860436772758&#93;
 &#91;2.9677861949066675, 0.15205611935372762, 25.600040161309696&#93;
 &#91;2.4046401797960795, 0.2914663505185634, 24.24373008707723&#93;
 &#91;1.9820054139405763, 0.46628657468365653, 22.964748583050085&#93;
 &#91;1.6788616460891923, 0.6565587545689172, 21.758445642263496&#93;
 &#91;1.4744010677851374, 0.8530017039412324, 20.62004063423844&#93;
</pre>



<pre class='hljl'>
<span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save_push</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
34.700 μs &#40;1010 allocations: 125.78 KiB&#41;
1000-element Array&#123;Array&#123;Float64,1&#125;,1&#125;:
 &#91;1.0, 0.0, 0.0&#93;
 &#91;0.8, 0.56, 0.0&#93;
 &#91;0.752, 0.9968000000000001, 0.008960000000000001&#93;
 &#91;0.80096, 1.3978492416000001, 0.023474005333333336&#93;
 &#91;0.92033784832, 1.8180538219817644, 0.04461448495326095&#93;
 &#91;1.099881043052353, 2.296260732619613, 0.07569952060880669&#93;
 &#91;1.339156980965805, 2.864603692722823, 0.12217448583728006&#93;
 &#91;1.6442463233172087, 3.5539673118971193, 0.19238159391549564&#93;
 &#91;2.026190521033191, 4.397339452147425, 0.2989931959555302&#93;
 &#91;2.5004203072560376, 5.431943011293093, 0.4612438424853632&#93;
 ⋮
 &#91;6.8089180814322185, 0.8987564841782779, 31.6759436385101&#93;
 &#91;5.6268857619814305, 0.3801973723631693, 30.108951163308078&#93;
 &#91;4.577548084057778, 0.13525687944525802, 28.545926978224173&#93;
 &#91;3.6890898431352737, 0.08257160199224252, 27.035860436772758&#93;
 &#91;2.9677861949066675, 0.15205611935372762, 25.600040161309696&#93;
 &#91;2.4046401797960795, 0.2914663505185634, 24.24373008707723&#93;
 &#91;1.9820054139405763, 0.46628657468365653, 22.964748583050085&#93;
 &#91;1.6788616460891923, 0.6565587545689172, 21.758445642263496&#93;
 &#91;1.4744010677851374, 0.8530017039412324, 20.62004063423844&#93;
</pre>


<p>This is because growth costs are amortized, meaning that when pushing, the size isn&#39;t increasing by one each time, but rather it&#39;s doing something like doubling, so that it&#39;s averaging O&#40;1&#41; cost to keep growing &#40;in theory the best is to do Golden ratio resizing, and long discussions can be had on this topic&#41;.</p>
<p>We can also look at what happens if we use matrices:</p>


<pre class='hljl'>
<span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save_matrix</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Matrix</span><span class='hljl-p'>{</span><span class='hljl-nf'>eltype</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>)}(</span><span class='hljl-n'>undef</span><span class='hljl-p'>,</span><span class='hljl-nf'>length</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>),</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-oB'>:</span><span class='hljl-p'>,</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u0</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-n'>n</span><span class='hljl-oB'>-</span><span class='hljl-ni'>1</span><span class='hljl-t'>
    </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-oB'>:</span><span class='hljl-p'>,</span><span class='hljl-n'>i</span><span class='hljl-oB'>+</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-oB'>:</span><span class='hljl-p'>,</span><span class='hljl-n'>i</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save_matrix</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
80.500 μs &#40;2001 allocations: 242.16 KiB&#41;
3×1000 Array&#123;Float64,2&#125;:
 1.0  0.8   0.752    0.80096   0.920338   …   1.98201    1.67886    1.4744
 0.0  0.56  0.9968   1.39785   1.81805        0.466287   0.656559   0.85300
2
 0.0  0.0   0.00896  0.023474  0.0446145     22.9647    21.7584    20.62
</pre>


<p>Where is this cost coming from? A large portion of the cost is due to the slicing on the <code>u</code>, which we can fix with a <code>view</code>:</p>


<pre class='hljl'>
<span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save_matrix_view</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Matrix</span><span class='hljl-p'>{</span><span class='hljl-nf'>eltype</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>)}(</span><span class='hljl-n'>undef</span><span class='hljl-p'>,</span><span class='hljl-nf'>length</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>),</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-oB'>:</span><span class='hljl-p'>,</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u0</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-n'>n</span><span class='hljl-oB'>-</span><span class='hljl-ni'>1</span><span class='hljl-t'>
    </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-oB'>:</span><span class='hljl-p'>,</span><span class='hljl-n'>i</span><span class='hljl-oB'>+</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-nd'>@view</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-oB'>:</span><span class='hljl-p'>,</span><span class='hljl-n'>i</span><span class='hljl-p'>]),</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save_matrix_view</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
36.801 μs &#40;1002 allocations: 132.89 KiB&#41;
3×1000 Array&#123;Float64,2&#125;:
 1.0  0.8   0.752    0.80096   0.920338   …   1.98201    1.67886    1.4744
 0.0  0.56  0.9968   1.39785   1.81805        0.466287   0.656559   0.85300
2
 0.0  0.0   0.00896  0.023474  0.0446145     22.9647    21.7584    20.62
</pre>


<p>Since we are only ever using single columns as a unit, notice that there isn&#39;t any benefit to keeping the whole thing contiguous, and in fact there are some downsides &#40;cache is harder to optimize because the longer cache lines are unnecessary, the views need to be used&#41;. Also, growing the matrix adaptively is not a very good idea since every growth requires both allocating memory and copying over the old values:</p>


<pre class='hljl'>
<span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save_matrix_resize</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Matrix</span><span class='hljl-p'>{</span><span class='hljl-nf'>eltype</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>)}(</span><span class='hljl-n'>undef</span><span class='hljl-p'>,</span><span class='hljl-nf'>length</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>),</span><span class='hljl-ni'>1</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-oB'>:</span><span class='hljl-p'>,</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u0</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-n'>n</span><span class='hljl-oB'>-</span><span class='hljl-ni'>1</span><span class='hljl-t'>
    </span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>hcat</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-nd'>@view</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-oB'>:</span><span class='hljl-p'>,</span><span class='hljl-n'>i</span><span class='hljl-p'>]),</span><span class='hljl-n'>p</span><span class='hljl-p'>))</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save_matrix_resize</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
1.860 ms &#40;2318 allocations: 11.69 MiB&#41;
3×1000 Array&#123;Float64,2&#125;:
 1.0  0.8   0.752    0.80096   0.920338   …   1.98201    1.67886    1.4744
 0.0  0.56  0.9968   1.39785   1.81805        0.466287   0.656559   0.85300
2
 0.0  0.0   0.00896  0.023474  0.0446145     22.9647    21.7584    20.62
</pre>


<p>So for now let&#39;s go back to the Vector of Arrays approach. One way to reduce the number of allocations is to require that the user provides an in-place non-allocating function. For example:</p>


<pre class='hljl'>
<span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>lorenz</span><span class='hljl-p'>(</span><span class='hljl-n'>du</span><span class='hljl-p'>,</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>α</span><span class='hljl-p'>,</span><span class='hljl-n'>σ</span><span class='hljl-p'>,</span><span class='hljl-n'>ρ</span><span class='hljl-p'>,</span><span class='hljl-n'>β</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>p</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-n'>α</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>σ</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-oB'>-</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]))</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-n'>α</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>ρ</span><span class='hljl-oB'>-</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>])</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>])</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-n'>α</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>β</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>])</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-n'>p</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-nfB'>0.02</span><span class='hljl-p'>,</span><span class='hljl-nfB'>10.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>28.0</span><span class='hljl-p'>,</span><span class='hljl-ni'>8</span><span class='hljl-oB'>/</span><span class='hljl-ni'>3</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Vector</span><span class='hljl-p'>{</span><span class='hljl-nf'>typeof</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>)}(</span><span class='hljl-n'>undef</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>similar</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u0</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-n'>n</span><span class='hljl-oB'>-</span><span class='hljl-ni'>1</span><span class='hljl-t'>
    </span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-n'>du</span><span class='hljl-p'>,</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
    </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-oB'>+</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>du</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-nf'>solve_system_save</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
1000-element Array&#123;Array&#123;Float64,1&#125;,1&#125;:
 &#91;1.0, 0.0, 0.0&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 ⋮
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
 &#91;-6.762970638284242, -11.790601841925698, 16.051639908100938&#93;
</pre>


<p>Oh no, all of the outputs are the same&#33; What happened? The problem is in the line <code>u&#91;i&#43;1&#93; &#61; du</code>. What we had done is set the save vector to the same pointer as <code>du</code>, effectively linking all of the pointers. The moral of the story is, you cannot get around allocating for all of these outputs if you&#39;re going to give the user all of the outputs&#33; It&#39;s impossible to not make all of these arrays, so if this is the case then you&#39;d have to:</p>


<pre class='hljl'>
<span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save_copy</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Vector</span><span class='hljl-p'>{</span><span class='hljl-nf'>typeof</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>)}(</span><span class='hljl-n'>undef</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>similar</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u0</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-n'>n</span><span class='hljl-oB'>-</span><span class='hljl-ni'>1</span><span class='hljl-t'>
    </span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-n'>du</span><span class='hljl-p'>,</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
    </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-oB'>+</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>copy</span><span class='hljl-p'>(</span><span class='hljl-n'>du</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-nf'>solve_system_save_copy</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
1000-element Array&#123;Array&#123;Float64,1&#125;,1&#125;:
 &#91;1.0, 0.0, 0.0&#93;
 &#91;0.8, 0.56, 0.0&#93;
 &#91;0.752, 0.9968000000000001, 0.008960000000000001&#93;
 &#91;0.80096, 1.3978492416000001, 0.023474005333333336&#93;
 &#91;0.92033784832, 1.8180538219817644, 0.04461448495326095&#93;
 &#91;1.099881043052353, 2.296260732619613, 0.07569952060880669&#93;
 &#91;1.339156980965805, 2.864603692722823, 0.12217448583728006&#93;
 &#91;1.6442463233172087, 3.5539673118971193, 0.19238159391549564&#93;
 &#91;2.026190521033191, 4.397339452147425, 0.2989931959555302&#93;
 &#91;2.5004203072560376, 5.431943011293093, 0.4612438424853632&#93;
 ⋮
 &#91;6.8089180814322185, 0.8987564841782779, 31.6759436385101&#93;
 &#91;5.6268857619814305, 0.3801973723631693, 30.108951163308078&#93;
 &#91;4.577548084057778, 0.13525687944525802, 28.545926978224173&#93;
 &#91;3.6890898431352737, 0.08257160199224252, 27.035860436772758&#93;
 &#91;2.9677861949066675, 0.15205611935372762, 25.600040161309696&#93;
 &#91;2.4046401797960795, 0.2914663505185634, 24.24373008707723&#93;
 &#91;1.9820054139405763, 0.46628657468365653, 22.964748583050085&#93;
 &#91;1.6788616460891923, 0.6565587545689172, 21.758445642263496&#93;
 &#91;1.4744010677851374, 0.8530017039412324, 20.62004063423844&#93;
</pre>


<p>which nullifies the advantage of the non-allocating approach. However, if only the end point is necessary, then the reduced allocation approach is helpful:</p>


<pre class='hljl'>
<span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_mutate</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-cs'># create work buffers</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>similar</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>);</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>copy</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-cs'># non-allocating loop</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-n'>n</span><span class='hljl-oB'>-</span><span class='hljl-ni'>1</span><span class='hljl-t'>
    </span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-n'>du</span><span class='hljl-p'>,</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
    </span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>du</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>du</span><span class='hljl-p'>,</span><span class='hljl-n'>u</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-nf'>solve_system_mutate</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
3-element Array&#123;Float64,1&#125;:
  1.4744010677851374
  0.8530017039412324
 20.62004063423844
</pre>


<p>Here we see a little trick: the line <code>u,du &#61; du,u</code> is swapping the pointer of <code>u</code> with the pointer of <code>du</code> &#40;since the value of the array is its reference&#41;. An alternative way to write the loop is:</p>



<pre class='hljl'>
<span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-n'>n</span><span class='hljl-oB'>-</span><span class='hljl-ni'>1</span><span class='hljl-t'>
  </span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-n'>du</span><span class='hljl-p'>,</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>.=</span><span class='hljl-t'> </span><span class='hljl-n'>du</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span>
</pre>


<p>which would compute <code>f</code> and then take the values of <code>du</code> and update <code>u</code> with them, but that&#39;s 3 extra operations than required, whereas <code>u,du &#61; du,u</code> will change <code>u</code> to be a pointer to the updated memory and now <code>du</code> is an &quot;empty&quot; cache array that we can refill &#40;this decreases the computational cost by ~33&#37;&#41;. Let&#39;s see what the cost is with this newest version:</p>


<pre class='hljl'>
<span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
38.901 μs &#40;1000 allocations: 109.38 KiB&#41;
3-element Array&#123;Float64,1&#125;:
  1.4744010677851374
  0.8530017039412324
 20.62004063423844
</pre>



<pre class='hljl'>
<span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_mutate</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
6.500 μs &#40;3 allocations: 336 bytes&#41;
3-element Array&#123;Float64,1&#125;:
  1.4744010677851374
  0.8530017039412324
 20.62004063423844
</pre>


<p>One last change that we could do is make use of StaticArrays. To do this, we need to go back to non-mutating, like:</p>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>StaticArrays</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>lorenz</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>α</span><span class='hljl-p'>,</span><span class='hljl-n'>σ</span><span class='hljl-p'>,</span><span class='hljl-n'>ρ</span><span class='hljl-p'>,</span><span class='hljl-n'>β</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>p</span><span class='hljl-t'>
  </span><span class='hljl-n'>du1</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-n'>α</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>σ</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-oB'>-</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]))</span><span class='hljl-t'>
  </span><span class='hljl-n'>du2</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-n'>α</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>ρ</span><span class='hljl-oB'>-</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>])</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>])</span><span class='hljl-t'>
  </span><span class='hljl-n'>du3</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-n'>α</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>β</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>])</span><span class='hljl-t'>
  </span><span class='hljl-nd'>@SVector</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-n'>du1</span><span class='hljl-p'>,</span><span class='hljl-n'>du2</span><span class='hljl-p'>,</span><span class='hljl-n'>du3</span><span class='hljl-p'>]</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-n'>p</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-nfB'>0.02</span><span class='hljl-p'>,</span><span class='hljl-nfB'>10.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>28.0</span><span class='hljl-p'>,</span><span class='hljl-ni'>8</span><span class='hljl-oB'>/</span><span class='hljl-ni'>3</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Vector</span><span class='hljl-p'>{</span><span class='hljl-nf'>typeof</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>)}(</span><span class='hljl-n'>undef</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u0</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-n'>n</span><span class='hljl-oB'>-</span><span class='hljl-ni'>1</span><span class='hljl-t'>
    </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-oB'>+</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-nf'>solve_system_save</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,</span><span class='hljl-nd'>@SVector</span><span class='hljl-p'>[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
1000-element Array&#123;SArray&#123;Tuple&#123;3&#125;,Float64,1,3&#125;,1&#125;:
 &#91;1.0, 0.0, 0.0&#93;
 &#91;0.8, 0.56, 0.0&#93;
 &#91;0.752, 0.9968000000000001, 0.008960000000000001&#93;
 &#91;0.80096, 1.3978492416000001, 0.023474005333333336&#93;
 &#91;0.92033784832, 1.8180538219817644, 0.04461448495326095&#93;
 &#91;1.099881043052353, 2.296260732619613, 0.07569952060880669&#93;
 &#91;1.339156980965805, 2.864603692722823, 0.12217448583728006&#93;
 &#91;1.6442463233172087, 3.5539673118971193, 0.19238159391549564&#93;
 &#91;2.026190521033191, 4.397339452147425, 0.2989931959555302&#93;
 &#91;2.5004203072560376, 5.431943011293093, 0.4612438424853632&#93;
 ⋮
 &#91;6.8089180814322185, 0.8987564841782779, 31.6759436385101&#93;
 &#91;5.6268857619814305, 0.3801973723631693, 30.108951163308078&#93;
 &#91;4.577548084057778, 0.13525687944525802, 28.545926978224173&#93;
 &#91;3.6890898431352737, 0.08257160199224252, 27.035860436772758&#93;
 &#91;2.9677861949066675, 0.15205611935372762, 25.600040161309696&#93;
 &#91;2.4046401797960795, 0.2914663505185634, 24.24373008707723&#93;
 &#91;1.9820054139405763, 0.46628657468365653, 22.964748583050085&#93;
 &#91;1.6788616460891923, 0.6565587545689172, 21.758445642263496&#93;
 &#91;1.4744010677851374, 0.8530017039412324, 20.62004063423844&#93;
</pre>



<pre class='hljl'>
<span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,</span><span class='hljl-nd'>@SVector</span><span class='hljl-p'>[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
7.934 μs &#40;2 allocations: 23.52 KiB&#41;
1000-element Array&#123;SArray&#123;Tuple&#123;3&#125;,Float64,1,3&#125;,1&#125;:
 &#91;1.0, 0.0, 0.0&#93;
 &#91;0.8, 0.56, 0.0&#93;
 &#91;0.752, 0.9968000000000001, 0.008960000000000001&#93;
 &#91;0.80096, 1.3978492416000001, 0.023474005333333336&#93;
 &#91;0.92033784832, 1.8180538219817644, 0.04461448495326095&#93;
 &#91;1.099881043052353, 2.296260732619613, 0.07569952060880669&#93;
 &#91;1.339156980965805, 2.864603692722823, 0.12217448583728006&#93;
 &#91;1.6442463233172087, 3.5539673118971193, 0.19238159391549564&#93;
 &#91;2.026190521033191, 4.397339452147425, 0.2989931959555302&#93;
 &#91;2.5004203072560376, 5.431943011293093, 0.4612438424853632&#93;
 ⋮
 &#91;6.8089180814322185, 0.8987564841782779, 31.6759436385101&#93;
 &#91;5.6268857619814305, 0.3801973723631693, 30.108951163308078&#93;
 &#91;4.577548084057778, 0.13525687944525802, 28.545926978224173&#93;
 &#91;3.6890898431352737, 0.08257160199224252, 27.035860436772758&#93;
 &#91;2.9677861949066675, 0.15205611935372762, 25.600040161309696&#93;
 &#91;2.4046401797960795, 0.2914663505185634, 24.24373008707723&#93;
 &#91;1.9820054139405763, 0.46628657468365653, 22.964748583050085&#93;
 &#91;1.6788616460891923, 0.6565587545689172, 21.758445642263496&#93;
 &#91;1.4744010677851374, 0.8530017039412324, 20.62004063423844&#93;
</pre>


<p>This is utilizing a lot more optimizations, like SIMD, automatically, which is helpful. Let&#39;s also remove the bounds checks:</p>


<pre class='hljl'>
<span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>lorenz</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>α</span><span class='hljl-p'>,</span><span class='hljl-n'>σ</span><span class='hljl-p'>,</span><span class='hljl-n'>ρ</span><span class='hljl-p'>,</span><span class='hljl-n'>β</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>p</span><span class='hljl-t'>
  </span><span class='hljl-nd'>@inbounds</span><span class='hljl-t'> </span><span class='hljl-k'>begin</span><span class='hljl-t'>
    </span><span class='hljl-n'>du1</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-n'>α</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>σ</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-oB'>-</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]))</span><span class='hljl-t'>
    </span><span class='hljl-n'>du2</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-n'>α</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>ρ</span><span class='hljl-oB'>-</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>])</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>])</span><span class='hljl-t'>
    </span><span class='hljl-n'>du3</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-n'>α</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>β</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>])</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-nd'>@SVector</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-n'>du1</span><span class='hljl-p'>,</span><span class='hljl-n'>du2</span><span class='hljl-p'>,</span><span class='hljl-n'>du3</span><span class='hljl-p'>]</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Vector</span><span class='hljl-p'>{</span><span class='hljl-nf'>typeof</span><span class='hljl-p'>(</span><span class='hljl-n'>u0</span><span class='hljl-p'>)}(</span><span class='hljl-n'>undef</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-nd'>@inbounds</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u0</span><span class='hljl-t'>
  </span><span class='hljl-nd'>@inbounds</span><span class='hljl-t'> </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-n'>n</span><span class='hljl-oB'>-</span><span class='hljl-ni'>1</span><span class='hljl-t'>
    </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-oB'>+</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-nf'>solve_system_save</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,</span><span class='hljl-nd'>@SVector</span><span class='hljl-p'>[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
1000-element Array&#123;SArray&#123;Tuple&#123;3&#125;,Float64,1,3&#125;,1&#125;:
 &#91;1.0, 0.0, 0.0&#93;
 &#91;0.8, 0.56, 0.0&#93;
 &#91;0.752, 0.9968000000000001, 0.008960000000000001&#93;
 &#91;0.80096, 1.3978492416000001, 0.023474005333333336&#93;
 &#91;0.92033784832, 1.8180538219817644, 0.04461448495326095&#93;
 &#91;1.099881043052353, 2.296260732619613, 0.07569952060880669&#93;
 &#91;1.339156980965805, 2.864603692722823, 0.12217448583728006&#93;
 &#91;1.6442463233172087, 3.5539673118971193, 0.19238159391549564&#93;
 &#91;2.026190521033191, 4.397339452147425, 0.2989931959555302&#93;
 &#91;2.5004203072560376, 5.431943011293093, 0.4612438424853632&#93;
 ⋮
 &#91;6.8089180814322185, 0.8987564841782779, 31.6759436385101&#93;
 &#91;5.6268857619814305, 0.3801973723631693, 30.108951163308078&#93;
 &#91;4.577548084057778, 0.13525687944525802, 28.545926978224173&#93;
 &#91;3.6890898431352737, 0.08257160199224252, 27.035860436772758&#93;
 &#91;2.9677861949066675, 0.15205611935372762, 25.600040161309696&#93;
 &#91;2.4046401797960795, 0.2914663505185634, 24.24373008707723&#93;
 &#91;1.9820054139405763, 0.46628657468365653, 22.964748583050085&#93;
 &#91;1.6788616460891923, 0.6565587545689172, 21.758445642263496&#93;
 &#91;1.4744010677851374, 0.8530017039412324, 20.62004063423844&#93;
</pre>



<pre class='hljl'>
<span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,</span><span class='hljl-nd'>@SVector</span><span class='hljl-p'>[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
5.517 μs &#40;2 allocations: 23.52 KiB&#41;
1000-element Array&#123;SArray&#123;Tuple&#123;3&#125;,Float64,1,3&#125;,1&#125;:
 &#91;1.0, 0.0, 0.0&#93;
 &#91;0.8, 0.56, 0.0&#93;
 &#91;0.752, 0.9968000000000001, 0.008960000000000001&#93;
 &#91;0.80096, 1.3978492416000001, 0.023474005333333336&#93;
 &#91;0.92033784832, 1.8180538219817644, 0.04461448495326095&#93;
 &#91;1.099881043052353, 2.296260732619613, 0.07569952060880669&#93;
 &#91;1.339156980965805, 2.864603692722823, 0.12217448583728006&#93;
 &#91;1.6442463233172087, 3.5539673118971193, 0.19238159391549564&#93;
 &#91;2.026190521033191, 4.397339452147425, 0.2989931959555302&#93;
 &#91;2.5004203072560376, 5.431943011293093, 0.4612438424853632&#93;
 ⋮
 &#91;6.8089180814322185, 0.8987564841782779, 31.6759436385101&#93;
 &#91;5.6268857619814305, 0.3801973723631693, 30.108951163308078&#93;
 &#91;4.577548084057778, 0.13525687944525802, 28.545926978224173&#93;
 &#91;3.6890898431352737, 0.08257160199224252, 27.035860436772758&#93;
 &#91;2.9677861949066675, 0.15205611935372762, 25.600040161309696&#93;
 &#91;2.4046401797960795, 0.2914663505185634, 24.24373008707723&#93;
 &#91;1.9820054139405763, 0.46628657468365653, 22.964748583050085&#93;
 &#91;1.6788616460891923, 0.6565587545689172, 21.758445642263496&#93;
 &#91;1.4744010677851374, 0.8530017039412324, 20.62004063423844&#93;
</pre>


<p>And we can get down to non-allocating for the loop:</p>


<pre class='hljl'>
<span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,</span><span class='hljl-nd'>@SVector</span><span class='hljl-p'>([</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>]),</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
5.867 μs &#40;1 allocation: 32 bytes&#41;
3-element SArray&#123;Tuple&#123;3&#125;,Float64,1,3&#125; with indices SOneTo&#40;3&#41;:
  1.4744010677851374
  0.8530017039412324
 20.62004063423844
</pre>


<p>Notice that the single allocation is the output.</p>
<p>We can lastly make the saving version completely non-allocating if we hoist the allocation out to the higher level:</p>


<pre class='hljl'>
<span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save!</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>n</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-nd'>@inbounds</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u0</span><span class='hljl-t'>
  </span><span class='hljl-nd'>@inbounds</span><span class='hljl-t'> </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-nf'>length</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>)</span><span class='hljl-oB'>-</span><span class='hljl-ni'>1</span><span class='hljl-t'>
    </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-oB'>+</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-n'>u</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-n'>u</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Vector</span><span class='hljl-p'>{</span><span class='hljl-nf'>typeof</span><span class='hljl-p'>(</span><span class='hljl-nd'>@SVector</span><span class='hljl-p'>([</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>]))}(</span><span class='hljl-n'>undef</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>solve_system_save!</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,</span><span class='hljl-nd'>@SVector</span><span class='hljl-p'>([</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>]),</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
4.957 μs &#40;0 allocations: 0 bytes&#41;
1000-element Array&#123;SArray&#123;Tuple&#123;3&#125;,Float64,1,3&#125;,1&#125;:
 &#91;1.0, 0.0, 0.0&#93;
 &#91;0.8, 0.56, 0.0&#93;
 &#91;0.752, 0.9968000000000001, 0.008960000000000001&#93;
 &#91;0.80096, 1.3978492416000001, 0.023474005333333336&#93;
 &#91;0.92033784832, 1.8180538219817644, 0.04461448495326095&#93;
 &#91;1.099881043052353, 2.296260732619613, 0.07569952060880669&#93;
 &#91;1.339156980965805, 2.864603692722823, 0.12217448583728006&#93;
 &#91;1.6442463233172087, 3.5539673118971193, 0.19238159391549564&#93;
 &#91;2.026190521033191, 4.397339452147425, 0.2989931959555302&#93;
 &#91;2.5004203072560376, 5.431943011293093, 0.4612438424853632&#93;
 ⋮
 &#91;6.8089180814322185, 0.8987564841782779, 31.6759436385101&#93;
 &#91;5.6268857619814305, 0.3801973723631693, 30.108951163308078&#93;
 &#91;4.577548084057778, 0.13525687944525802, 28.545926978224173&#93;
 &#91;3.6890898431352737, 0.08257160199224252, 27.035860436772758&#93;
 &#91;2.9677861949066675, 0.15205611935372762, 25.600040161309696&#93;
 &#91;2.4046401797960795, 0.2914663505185634, 24.24373008707723&#93;
 &#91;1.9820054139405763, 0.46628657468365653, 22.964748583050085&#93;
 &#91;1.6788616460891923, 0.6565587545689172, 21.758445642263496&#93;
 &#91;1.4744010677851374, 0.8530017039412324, 20.62004063423844&#93;
</pre>


<p>It is important to note that this single allocation does not seem to effect the timing of the result in this case, when run serially. However, when parallelism or embedded applications get involved, this can be a significant effect.</p>
<h2>Discussion Questions</h2>
<ol>
<li><p>What are some ways to compute steady states? Periodic orbits?</p>
</li>
<li><p>When using the mutating algorithms, what are the data dependencies between different solves if they were to happen simultaniously?</p>
</li>
<li><p>We saw that there is a connection between delayed systems and multivariable systems. How deep does that go? Is every delayed system also a multivariable system and vice versa? Is this a useful idea to explore?</p>
</li>
</ol>


        <HR/>
        <div class="footer">
          <p>
            Published from <a href="dynamical_systems.jmd">dynamical_systems.jmd</a>
            using <a href="http://github.com/JunoLab/Weave.jl">Weave.jl</a> v0.10.3 on 2020-09-15.
          </p>
        </div>
      </div>
    </div>
  </div>
</BODY>

</HTML>
